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The Terminal Area Simulation System (TASS) is a large eddy scale atmospheric flow 
model with extensive turbulence and microphysics packages. It has been applied 
successfully in the past to a diverse set of problems ranging from prediction of severe 
convective events (Proctor et al. 2002), tracking storms and for simulating weapons effects 
such as the dispersion and fallout of fission debris (Bacon and Sarma 1991), etc. More 
recently, TASS has been used for predicting the transport and decay of wake vortices behind 
aircraft (Proctor 2009). An essential part of the TASS model is its comprehensive 
microphysics package, which relies on the accurate computation of microphysical scalar 
transport. This paper describes an evaluation of the Leonard scheme implemented in the 
TASS model for transporting microphysical scalars. The scheme is validated against 
benchmark cases with exact solutions and compared with two other schemes - a Monotone 
Upstream-centered Scheme for Conservation Laws (MUSCL)-type scheme after van Leer 
and LeVeque’s high-resolution wave propagation method. Finally, a comparison between 
the schemes is made against an incident of severe tornadic super-cell convection near Del 
City, Oklahoma. 


I. Introduction 

M icrophysics packages in atmospheric flow models depend on the accurate computation of scalar transport of 
water vapor, ice, snow, and hail, etc. These microphysical scalars are fully coupled to the dynamical solver in 
the sense that they provide heat sources/sinks which, in turn, alter the flow field accordingly. Any algorithm for 
solving the scalar transport equation should preferably have the following desirable properties: conservation, 
positivity, shape preservation in addition to accuracy and stability. Traditionally, the dynamical cores of 
atmospheric flow models have been based on centered schemes which can introduce false negatives in the solution 
of scalar transport equation. Therefore, for the advection of scalars, positive-definite schemes are used. A large 
number of such schemes have been suggested over the years and an excellent review is given by LeVeque (1996). 
Some recent positive-definite schemes include the flux-corrected transport FCT algorithm (Zalesak 1979), 
Leonard’s third-order upwind-biased scheme (Leonard 1988; Leonard et al. 1993; Leonard et al. 1995), the 
Multidimensional Positive Definite Advection and Transport Algorithm (MPDATA) suggested by Smolarkiewicz 
(1984), higher-order MUSCL-type schemes after van Leer (1979), high-resolution wave propagation methods 
(LeVeque 1996), and weighted area flux (WAF) and the weighted, essentially non-oscillatory (WENO) schemes 
(Titarev and Toro 2005: Toro 1999). 

In this study, the Leonard scheme is compared with a MUSCL-type scheme and the high-resolution wave 
propagation method suggested by LeVeque. Both schemes use limiters to enforce the total variation diminishing 
(TVD) condition to maintain positive-definiteness of the solution. The objective was to down select a positive 
definite scheme for the TASS model. TASS is a large eddy scale atmospheric flow model with comprehensive 
turbulence and microphysics schemes. It is a fully parallelized code (using MPI) that exhibits excellent scaling 
across processors. Computational efficiency of the scheme and the ease of implementation within the existing 
TASS code infrastructure were two other important factors in selecting the appropriate scheme in addition to 
desirable properties of the advection solver mentioned earlier. In the following sections, the three numerical 
schemes (Leonard scheme, MUSCL-type scheme and the high-resolution wave propagation method) are described 
in detail. The schemes are analyzed and evaluated against benchmark cases in two and three dimensions with exact 
solutions and against a real incident of severe tornadic super-cell convection near Del City, Oklahoma. 
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II. Numerical Schemes 

The two-dimensional unsteady scalar transport equation in conservative form and with the assumption of 
incompressibility for any microphysical scalar quantity q (vapor, ice, snow, hail, etc.) can be written as: 
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where, 


F = qu, G - qv . 
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In Eq. (l)-(2), F and G are the fluxes in the x and y directions, u and v are the two velocity components, k is the 
eddy diffusivity and S q is the source term. The focus of this study is on scalar transport and therefore the diffusion 
and source terms are neglected. Eq. (l)-(2) can be written in the discrete flux-difference form as: 
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Eq. (3) can be re-written as: 
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where, the update in time n + 1 of the conserved scalar, q.., is obtained by summing up all the fluxes in the x 

direction (F i+V2j and F._ V2j ) and in they direction ( G. /+1/2 and G ij _ V2 ) across the control surfaces of each control 

volume (see Figure 1). In the following sections, the Leonard scheme, a MUSCL-type scheme and LeVeque’s high- 
resolution wave propagation method are described in some detail. It should be noted that although the schemes are 
described with the help of a 2D stencil and most of the test cases are two dimensional, the final implementation of 
the algorithms in the TASS code is in three dimensions. 
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Figure 1. Control volume in a two-dimensional Cartesian mesh. TASS uses the Arakawa-C grid staggering. 
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A. Leonard Scheme 

The Leonard scheme is a higher-order, fully multidimensional, conservative, finite volume scheme (Leonard 
1988; Leonard et al. 1993; Leonard et al. 1995). It is formally third-order accurate in space and exhibits minimal 
numerical dispersion, diffusion and phase errors. Leonard’s scheme can be written as: 
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The flux F ._ V2j , across the left face of the control volume q tj (see Figure 2) is given by: 
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where, c = i/ — and c =v — are the Courant numbers in the x and y directions, respectively. Fluxes F . +V2 . , 
Ax J Ay ,J 

G ij _ l/2 , and G ij+ll2 can be constructed in a similar manner. 

The basic stencil for calculating the Leonard flux (see Figure 2) requires seven points. Versions of the algorithm 
that enforce the TVD condition have been suggested by Leonard (1988) but suffer from directional errors due to 
their unidirectional construction, as pointed out by Thubum (1996). Thuburn’s alternative is a multidimensional 
limiter which performs much better than the universal limiter suggested by Leonard. It is, however, difficult to 
implement within the existing parallelized code structure and also computationally expensive. In this study the 
dispersion errors were corrected by using a simple methodology - the negative values were set to zero followed by a 
re-normalization step to enforce mass conservation. The correction is applied every seventh time step. The results 
from this corrected version of the scheme are denoted by LeonardC in the results section. 
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Figure 2. Stencil for calculating the third-order upwind-biased Leonard flux. 
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B. MUSCL-type Scheme with Barth- Jesperson Slope Limiter 

The description of a Monotone Upstream-centered Scheme for Conservation Laws (MUSCL)-Type scheme is 
given in this section. Higher-order accuracy in space for the Godunov method (Godunov 1959) can be achieved by 
constructing piecewise linear data from cell averages (van Leer 1979): 


Hl - Hu 




( 7 ) 


H R =<lir+(Xfaee ~ X ir) L face + O face ~ T ir) L face > ( 8 ) 

where, q L is the extrapolated value of the conserved quantity q on the left side of the face of the control volume, the 
subscript il is used for the cell center quantities in the cell on the left of the face, e.g., qn is the cell-averaged value of 
the conserved quantity q and (xu,yu) are the cell center coordinates of the cell on the left. Similarly, q R is the 
extrapolated value of the conserved quantity q on the right side of the face, the subscript ir is used for the cell center 
quantities in the cell on the right of the face, e.g., q ir is the cell-averaged value of the conserved quantity q and 
(Xi n y ir ) are the cell center coordinates of the cell on the right. The point of intersection between the face and the line 
connecting the centers of the two cells on either side of the face is denoted by (xf ace ,yf ace ). Lf ace is the limiter on the 
gradient to ensure a monotonic solution. For higher spatial accuracy, these extrapolated values are used in the flux 
calculations instead of cell-centered averages. On a Cartesian mesh, implementation is straightforward, e.g., q it is 
defined by q (i-1 , j ) and q ir by q (i , j ) . See Figure 3. The gradients can be calculated using the 2 nd or 4 th order 
difference formulae. The flux F ._ ll2j , across the left face of the control volume q tj (Figure 3) is given by: 

= u i,flA y ( 9 ) 

The piecewise linear reconstruction of data is bounded by enforcing the following condition (Harten 1983: Barth 
and Jesperson 1989), 
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where, Ny are the neighbors of the cell (i, j) . The neighborhood cells can be defined in different ways. In this 
study, it was found that using cells (i-l,j), (i-1, j-1) , (i-1, j+1) , (i, j-1) , (i,j+l), (i+1, j) , 

(i+1, j-1) and (i+1, j+1) as the neighborhood around the cell (i, j) gave good results. The limiter Lf ace can 
now be determined by: 
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where, q.. is the cell-averaged value and cf ace is the extrapolated value on the face of the cell. Four values of 

Lf ace are obtained for each cell (one for each face) from Eq. (13) and the minimum of the four is used to limit the 
gradient at the cell center. The gradients are calculated using second-order differences. A two-stage explicit Runge- 
Kutta time integration scheme (Jameson et al. 1981) can be used to march the solution in time. Let, 


R 



Then the two-stage Runge-Kutta time-marching scheme can be written as: 
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Where, oti = 1/2, and a 2 = 1 are the Runge-Kutta constants. 
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Figure 3. Stencil for calculating the MUSCL flux in two dimensions. 

C. LeVeque’s High-resolution Wave Propagation Method 

In the LeVeque scheme (Leveque 1996; LeVeque 2002) first the Riemann problem is solved across each face to 
calculate the waves and their speeds. In the 2D scalar transport equation, there is only one wave that propagates 
across each control surface. For example, in the x direction, the wave W ._ l/2J , and its speed s._ U2J are given by: 
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The fluctuations A + Aq._ l/2 propagating across the interface can be computed as follows (in the x direction): 
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The correction fluxes are computed as follows: 
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A second Riemann solver in the transverse direction uses the fluctuations A ± Aq l _ U2J to find the up and down 
going fluctuations, e.g., A + Aq._ l/2J is used to calculate: 
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These are used to update the correction fluxes G t ,, 2 and G. . +1/2 respectively. Similarly, ^4 Ag, 
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The second-order correction is given by the following equation: 
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where, JT._ 1/2 / is the limited version of the wave, JT._ 1/2 . . A similar expression can be obtained for G._ ]/2/ and 

other fluxes across the control surfaces of each control volume. This results in a general and a very robust algorithm 
which can be used to solve any hyperbolic conservation law as long as its eigenstructure is known. The 
implementation of the algorithm however is involved (especially given the parallel constructs of the TASS model) 
and it is relatively expensive in computational terms. The coding of the scheme for this study was greatly aided by 
closely following the CLAWPACK (LeVeque 2003) software implementation. 
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III. Results 


The Leonard scheme is validated against different benchmark cases in two dimensions - rotating Gaussian 
function, advection of square pulse in a uniform diagonal flow, Doswell’s smooth frontogenesis and 
Smolarkiewicz’s deformational flow. In three dimensions, the scheme is validated against LeVeque’s deformational 
flow. Noye-Tan test for advection-diffusion is also performed to evaluate the diffusion flux discretization. The 
results from the Leonard scheme are also compared against solutions obtained for the first-order upwind scheme, the 
high-resolution wave propagation method suggested by LeVeque and a MUSCL-type scheme with slope limiter. 

Finally, the Leonard scheme in the TASS model is compared against the MUSCL implementation in the TASS 
model using the Del City case of severe convection using real sounding data. 

Cases A through F were run using a standalone code which was written to evaluate the different schemes. Cases 
G and H were run using the schemes implemented in the TASS model. The computational domains in Cases A 
through F have arbitrary spatial units in both the x- and y-direction. The following metrics were used in the error 
analysis: The RMS error is given by: 
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The error is defined in terms of the L r norm as: 
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where, ncells is the total number of cells in the mesh. The diffusion error was found by subtracting the computed 
tracer maxima from the exact value of the tracer maxima (Iselin et al., 2002): 

E diffusion = ma .x(<T“ ) - maxfo— ) . (32) 


The RMS and Lj errors were used in a mesh convergence study to determine the convergence rate of each 
scheme. The convergence rate, p of a numerical scheme is given by: 
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where, E 2h /E h is the ratio of error between meshes of resolutions h and 2h. 


A. Rotating Cone Test 

The rotating cone case was run for different schemes in a mesh convergence study to compare the convergence 
rates. The domain for this case was bounded within [0,100] x [0,100]. Three different meshes with a resolution of 
0.5, 1 and 2, respectively, were used in the convergence study. 

A smooth Gaussian cone function (radius, r b = 10) was initialized, centered at (x c , y c ) = (50, 0.75y max ): 
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The rotational flow field was defined as: 
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u(x,y) = -a)(y-y 0 ). 


(36) 


v(x,y)= co(x-x o ), (37) 

where, u{x,y) and v(x,y) are the velocities in the x and y directions, respectively. (x 0 , y 0 ) are the coordinates of the 
mesh center. The constant angular velocity co was set to 0.1. The simulation was run for 62.832s, which is the time 
taken by the cone to complete one revolution = In / co. The rotating cone tests were performed for four different 
schemes: the Leonard scheme; Leonard’s scheme with correction for negative values (LeonardC); a MUSCL-type 
scheme with the Barth-Jesperson slope limiter (MUSCL-WL); and LeVeque’s high-resolution wave propagation 
scheme. 

Farfield boundary conditions were defined with the help of ghost cells in each case. For Leonard and MUSCL 
schemes one extra layer of cells was used to define boundary conditions and two layers of cells for setting boundary 
conditions were used for the LeVeque scheme. For the MUSCL scheme, the solution was marched in time using a 
two-stage explicit Runge-Kutta scheme (Jameson et al., 1981). A CFL criterion of 0.25 was set for Leonard’s 
scheme and 0.9 for all other schemes. Different limiter functions were tested with the LeVeque scheme ( minmod , 
superbee , van Leer, van Albada and monotonized-centered). The monotonized-centered limiter gave the best 
performance and was used for computations shown in this paper. 

Figure 4 shows the initial conditions for this test case. The convergence rates for different schemes (p L \ and p m s ) 
are tabulated in Table 1. The Leonard and the LeonardC schemes have the fastest convergence rates both in terms 
of RMS error. In terms of the L r norm, the MUSCL scheme is comparable to the Leonard scheme. The LeVeque 
scheme had the slowest convergence rate when compared to the Leonard and the MUSCL-type scheme with limiter. 
The decrease in errors as a function of the mesh resolution is plotted in Figure 5. The timings for the different 
schemes as a function of error reduction and mesh resolution are also plotted in Figure 5. Since both the MUSCL- 
type and the LeVeque schemes run at a higher CFL number, their timing is also less than that of the Leonard 
scheme. For the sake of completeness, one additional run was made with the MUSCL-scheme with limiter at the 
same CFL number as the Leonard scheme (CFL = 0.25) and the results are plotted in Figure 5 alongside other runs. 
It can be seen that at the same CFL number, the Leonard scheme is comparable to the MUSCL scheme in terms of 
computational efficiency. LeVeque’s high-resolution wave propagation method was found to be computationally 
more expensive than the MUSCL-type scheme. The different schemes were timed using the Linux command time. 

The concentration contours of the scalar field after one complete revolution for the different schemes are shown 
in Figure 6 for a mesh resolution of 1. The concentration contours of the scalar field after one complete revolution 
for the different schemes are shown in Figure 7 for a mesh resolution of 0.5. The Leonard scheme and the Leonard 
scheme with correction for negative values (LeonardC) have better shape preserving property compared to the other 
two schemes for this particular test case. The correction for negative values in the Leonard scheme results in a 
scheme that is slightly more diffusive than the Leonard scheme but still has a higher maximum compared to the 
other two schemes. LeVeque’s high-resolution wave propagation method is more diffusive than the MUSCL-type 
scheme. 



Figure 4. Initial conditions: scalar [min = 0 : max = 0.987] (left panel), w-velocity [min = -4.95 : max = 4.95] 
(middle panel), and v-velocity [min = -4.95 : max = 4.95] (right panel). The axis in both horizontal and 
vertical have arbitrary units. 
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Table 1: Convergence Rates 
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Figure 5. Reduction in error for different schemes with increasing mesh resolution and associated 
computational cost. In the plots, dx is the grid spacing in arbitrary units. Time has units of seconds. 
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Figure 6. The results of a rotating Gaussian cone test using, Leonard scheme (top left), Leonard scheme with 
correction (top right), MUSCL with Barth-Jesperson slope limiter (bottom left), and LeVeque scheme with 
limiter (bottom right) after one revolution (Time = 62.832s). The grid spacing was set to dx = dy = 1 
(arbitrary units). 100 x 100 cells. The axis in both horizontal and vertical have arbitrary units. 
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Figure 7. The results of a rotating Gaussian cone test using, Leonard scheme (top left), Leonard scheme with 
correction (top right), MUSCL with Barth-Jesperson slope limiter (bottom left), and LeVeque scheme with 
limiter (bottom right) after one revolution (Time = 62.832s). The grid spacing was set to dx = dy = 0.5 
(arbitrary units). 200 x 200 cells. The axis in both horizontal and vertical have arbitrary units. 


11 

American Institute of Aeronautics and Astronautics 



B. Square Pulse in Diagonal Flow 

A square function of unit height was placed in a background scalar field of 0.1. A uniform diagonal flow field 
was imposed with u(x,y) = 0.5 and v(x,y ) = 1.0. The domain was bounded within [0,1] x [0,1] and a total of 100 
cells were used in both the x and y directions. Periodic boundary conditions were used in all directions. Figure 8 
shows the computed solutions at t = 2s. The results from the first-order upwind scheme are shown as a reference. 
Both the Leonard and LeVeque schemes have better shape preservation because both schemes take into account 
cross derivative terms. The MUSCL scheme with Barth-Jesperson limiter, on the other hand, maintains the minima 
and maxima throughout the simulation, without any undershoots or overshoots. In the LeVeque scheme, there is a 
small undershoot, whereas the Leonard scheme shows errors due to the presence of discontinuities. 



Figure 8. First-order upwind solution (top left), Leonard scheme without correction (top right), MUSCL with 
Barth-Jesperson slope limiter (bottom left), and LeVeque scheme with limiter (bottom right) after Time = 2s. 
The grid spacing was set to dx = dy = 0.005 (arbitrary units). 100 x 100 cells. Periodic boundary conditions 
were used in this calculation. The axis in both horizontal and vertical have arbitrary units. 
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C. Doswell’s Frontogenesis 

The simulation of DoswelTs frontogenesis problem (Doswell 1984; Holm 1995) is presented in this section. 
Doswell’s idealized model describes the interaction of a nondivergent vortex with an initially straight frontal zone. 
The domain in this case was bounded within [-4,4] x [-4,4]. The simulation was run for t = 4 seconds. The flow 
field for the Doswell test is defined as follows: 


u(x,y)= 7 /' , 

r /max 

(38) 

v(x,y)=--L 

r f 

J max 

(39) 

where, r is the distance from any given point to the origin of the coordinate system (i.e. the point: x = 0 andy = 0), 
f max = 0.385 is the maximum tangential velocity and f is given by: 

tanh(r) 
cosh 2 (r) 

The evolution of tracer field in time t , is given by the exact solution: 
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where, 5 is the characteristic width of the frontal zone. The value of S was set to 2. The initial scalar field can be 
obtained from Eq. (41) by setting t = 0. Figure 9 shows the initial conditions for the test case. 

The comparison of different schemes at t = 4s is shown in Figure 10 along with the exact solution. The 
computed results from different schemes are also compared with the exact solution along the mesh centerline in 
Figure 1 1 for different mesh resolutions. The computed results for all schemes are in excellent agreement with the 
exact solution (see, e.g ., Holm 1995, Ahmad et al. 2006). The behaviors of the MUSCL-type scheme and the high- 
resolution wave propagation method are quite similar. The Leonard scheme is less diffusive than both the MUSCL- 
type scheme and the high-resolution wave propagation method and is able to fully capture the distortion of the front, 
both inside and outside of the stationary vortex. 




Figure 9. DoswelPs Frontogenesis. Initial conditions: scalar [min =-0.962 : max = 0.962] (left panel), u- 
velocity [min = -0.995 : max = 0.995] (middle panel), and v-velocity [min = -0.995 : max = 0.995] (right panel). 
The axis in both horizontal and vertical have arbitrary units. 
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Figure 10. Doswell’s Frontogenesis. Comparison of different computed solutions with the exact solution 
for a smooth frontogenesis at time = 4s. The grid spacing was set to dx = dy = 1 (arbitrary units). The axis in 
both horizontal and vertical have arbitrary units. 
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Figure 11. Doswell’s Frontogenesis. Comparison of different computed solutions with the exact solution 
for a smooth frontogenesis at time = 4s. Centerline data between x = -4 and x = 4 at y = 0. The first-order 
upwind solution is shown as a reference in the top panel. The axis in both horizontal and vertical have 
arbitrary units. 


15 

American Institute of Aeronautics and Astronautics 


D. Smolarkiewicz’ s Deformational Flow 

Smolarkiewicz’s deformational flow (Smolarkiewicz 1982; Staniforth et al. 1987) is described in this section. 
The domain was bounded within [0,100] x [0,100]. A scalar cone with a height of 1 unit and radius, R = 15 units 
was initialized in the middle of the domain: 


q(x,y) = \-^- if r <R , 
K 


(43) 


where, 


r = Tjix-X'Y+iy-yJ 2 , (44) 

and, (x 0 , y 0 ) is the center of the mesh. The flow field for the deformation test consists of sets of symmetrical vortices 
and is given by: 


u{x , y) - Ak sin {kx) sin(Ay) ; v(x, y) - Ak cos {kx) cos {ky ) , (45) 

where, u(x,y) and v(x,y) are again the velocities in the x and y directions, respectively, k = 47r/L , 
A = 8 and L = 100 units. The initial conditions for this test case are shown in Figure 12. 

Figure 13 shows the scalar distribution at t = T/50 for different schemes, where, T = 2637.6 seconds is the final 
time of integration used in Smolarkiewicz (1982). The domain maximum and minimum value of the scalar field is 
also shown in the figure. The Leonard scheme with correction for negative values is again denoted by LeonardC. 
All schemes produce a symmetric result and except for the MUSCL-type scheme, negative values appear in the 
solution. 

Figures 14 and 15 show the comparison with Staniforth’s analytical solution for scalar values between x = 20 
and 80, for y = 50. The Staniforth solution is computed numerically and requires an input of sampling interval. For 
the comparison shown in Figures 14-15, a sampling interval of 0.1 was used for the analytical solution. Staniforth et 
al. (1987) have discussed this test case in detail. They point out that for a mesh resolution of 1 used in 
Smolarkiewicz (1982), the numerical solution is valid only for time < T/50. After time > T/50 the features of the 
tracer field become too small to be effectively captured by a mesh resolution of 1, normally used for this test. Even 
with higher resolution, some of the features (the spike in the middle) are sub-grid scale for the schemes to resolve. 
Overall, the computed results are in excellent agreement with the analytical solution (Staniforth 1987) and 
previously reported results of the benchmark (e.g., Ahmad et al. 2006). The Leonard and LeonardC schemes are 
both less diffusive than the MUSCL-type scheme and the high-resolution wave propagation method. 



20 40 60 80 20 40 60 80 20 40 60 80 


Figure 12. Smolarkiewicz Flow. Initial conditions: scalar [min = 0 : max = 0.976] (left panel), w-velocity 
[min = -1 : max = 1] (middle panel), and v-velocity [min = -1 : max = 1] (right panel). The axis in both 
horizontal and vertical have arbitrary units. 


16 

American Institute of Aeronautics and Astronautics 





















Figure 13. Smolarkiewicz Flow. The scalar field at time = T/50s into the simulation is shown for different 
schemes. T = 2637.6s. The grid spacing was set to dx = dy = 1 (arbitrary units). The axis in both horizontal 
and vertical have arbitrary units. 
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Figure 14. Smolarkiewicz Flow. Comparison of the scalar field at time = T/50s into the simulation with 
exact solution. Centerline data is shown for x = 20 to x = 80 at y = 50. T = 2637.6s. In the top left panel, the 
results from the first-order upwind scheme are also shown for reference. The axis in both horizontal and 
vertical have arbitrary units. 
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E. Noye-Tan Test 

The focus of this study is on the evaluation of different advection algorithms. The final implementation in the 
TASS model however requires the computation of the advection-diffusion equation for various microphysical 
scalars. For the sake of completeness, the results of the test case for the advection-diffusion equation proposed by 
Noye and Tan (1989) are presented in this section. The test case describes the diffusion of an initial Gaussian pulse 
as it is advected along a straight line. The domain was bounded within [0,3] x [0,3] and 100 cells were used in each 
direction. The analytical solution of the unsteady advection-diffusion for the Noye-Tan test case is given by: 


q(x,y,t) 


1 ( x-ut-x 0 ) 2 ( 7 -V/- 7 J 2 

4t + \ p |_ k r (4t + l) k y {At + 1) 


(46) 


where, k x = k y = 0.01 are the diffusion coefficients and u = v = 0.8 are the velocities in the x and y directions, 
respectively. x 0 = y 0 = 0.5 is the center of the initial tracer distribution. The initial conditions and the Dirichlet 
boundary conditions can be obtained from the analytical solution by setting t = 0. In the case of LeonardC scheme, 
the third-order diffusion flux suggested by Leonard (Leonard et al. 1993) was used. Second-order differences were 
used to add the diffusion term in the MUSCL-type scheme. After the advection step, the scalar is updated by adding 
the contribution due to diffusion in both the x and y directions: 


At 


-c)-Kh M -C'J] 


(47) 


In Eq. (47), £ 2 , & 2 , etc. are the eddy diffusivities at control surfaces. The Ax 2 and the Ay 2 terms in Eq. 

(47) put a severe time step restriction on the diffusion step. The simulation was run for t = 2s. The initial conditions 
and the computed solution using LeonardC scheme at t = 2s are shown in Figure 16. A comparison between the 
exact solution and the computed solution along the mesh diagonal is shown in Figure 17 for the MUSCL-type and 
LeonardC schemes at different times. 



Figure 15. Noye-Tan Test, scalar initial profile (left panel), solution at Time = 2 using LeonardC scheme 
(right panel). The axis in both horizontal and vertical have arbitrary units. 


19 

American Institute of Aeronautics and Astronautics 





cd 

o 

CO 



d 


Figure 16. Noy e-Tan Test. Comparison of the LeonardC scheme with the MUSCL-type scheme. The 

computed solution and the exact solution are shown along the mesh diagonal at different times: Time = 0 
(green), 1 (red) and 2 (blue). Exact solution is given by solid lines and the computed solutions are represented 
by points. The axis in both horizontal and vertical have arbitrary units. 
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F. LeVeque’s Three Dimensional Deformation Flow 

In this section, the Leonard scheme is compared with the MUSCL-type scheme in three dimensions using a test 
case suggested by LeVeque (1996). This is an interesting and relatively severe test case - the flow field is truly 
three-dimensional as well as time-dependent and, in addition, the scalar field is characterized by a discontinuous 
function. The domain for this test case was bounded within [0,1] x [0,1] x [0,1] and 40 cells were used in each 
direction. The time-dependent flow field is given by: 


u(x,y,z ) = 2sin 2 (;rx)sin(2;ry)sin(2;rz)g(£) 

(48) 

v(x,y,z) = -sin(2;rx)sin 2 (;ry)sin(;rz)g(f) 

(49) 

w(x,y,z ) = -sin(2;rx)sin(;ry)sin 2 (;rz)g(f) 

(50) 


The time dependence, g(£) of the flow field is defined as: 


g(t) = COS 



(51) 


The constant T in Eq. (51) is set to 1.5. The flow reverses direction at t = 772, therefore the initial scalar field 
should be recovered at t = T. The initial scalar field is set according to the following relation: 


q(x,y,z) 


f 1, if x < 1 / 2, 
[0, x > 1/2. 


(52) 


Figure 18 shows the initial three dimensional velocity field at in the x-y plane at z = 0.4125. The Leonardo 
scheme and the MUSCL-type scheme with slope limiter are compared in Figure 19. The cross section at in the x-y 
plane at z = 0.4125 is shown at t= T/2 and t = T for the two schemes. Since, there is a correction for negative values 
in the LeonardC scheme the minima is maintained at 0, the maximum however overshoots the original value by 
0.114 at t = T/2 and by 0.057 at simulation end time. The oscillation in the maximum value is to be expected 
because the Leonard scheme is not TVD. The MUSCL-type scheme on the other hand maintains the sharp interface 
without any overshoots or undershoots throughout the simulation. 



Figure 17. LeVeque’s Deformation Flow. Initial conditions: w-velocity [min = -1.04 : max = 1.04] (left 
panel), v-velocity [min = -0.519 : max = 0.519] (center panel), and w-velocity [min = -0.919 : max = 0.919] 
(right panel). Cross section in the x-y plane at z = 0.4125 is shown in the figure. The axis in both horizontal 
and vertical have arbitrary units. 
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Figure 18. LeVeque’s Deformation Flow. Comparison of the LeonardC scheme with the MUSCL-type 
scheme. LeonardC at time = 172 (top left); MUSCL at time = 172 (top right); LeonardC at time = T (bottom 
left) and MUSCL at time = T (bottom right), x-y plane is shown at z = 0.4125. The final time of the 
simulation is, T = 1.5s. The axis in both horizontal and vertical have arbitrary units. 
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G. Idealized Microburst 

A simulation of an idealized microburst was run to test the implementation of MUSCL-type scheme with Barth- 
Jesperson limiter in the TASS model. The domain was bounded within [0,6km] x [0,6km] x [0,3km]. The mesh 
resolution in all three directions was set to 60m and the simulation was run for a final time of 400s. TASS model 
was initialized by a single sounding which was horizontally homogeneous (Figure 19) and run with full 
microphysics. The initial velocity field was set to zero. See Proctor (1988) for further details on this test case. The 
Orlanski symmetry ramping function is introduced at the model domain top to specify precipitation: 


, , * , 17/150 

q h {x,y,t) = q m y- e(r xj j 


t < 1505 
1505 <t < 4005 


(53) 


In Eq. (53), q h controls the maximum amount of precipitation as a function of space and time that is introduced 
at the top of the computational domain. q max = 6g/m 3 is the maximum amount of precipitation, 

r - yj(x-x o ) 2 +(y~y o ) 2 is the distance from the center of precipitation initialization (x 0 , yo), and cr is the sigma 

value of precipitation distribution at the top of the domain. A comparison of the Leonard scheme with the MUSCL- 
type for this test case is given in Figures (20)-(22). Figures 20 and 21 show the comparison of hail and rain fields 
after 200s and 400s into the simulation, respectively. The x-z plane at y = y max /2 is displayed. Overall, the 
differences between the two simulations are small. The MUSCL-type scheme maintains the positive-definiteness of 
the scalar fields whereas small negative values appear in the Leonard solution. The Leonard scheme is also slightly 
more diffusive than the MUSCL-type scheme. Figure 22 shows the comparison between the two schemes for the 
vertical velocity and hail at time = 400s into the simulation. A zoomed-in view of the x-y plane is shown at z = 
z max /2. Again the main difference between the two schemes is in the preservation of positivity. The MUSCL-type 
scheme is able to maintain the positive definiteness of the scalars. These differences in scalars also alter the flow 
field slightly. The maximum and minimum values of vertical velocities at the same x-y plane shown in Figure 22 
are different for the Leonard and MUSCL-type scheme. Similar differences in the horizontal velocities were also 
observed. In terms of computational efficiency the MUSCL-type scheme was slightly faster than the Leonard 
scheme. 
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Figure 19. Idealized Microburst: Initialization sounding (left). Initial phase of the flow development and 
formation of hail at time = 20s into the simulation is shown in the right panels. Vertical velocity in m/s (top 
right) and hail in kg/m 3 (bottom right), x-z plane is shown at y = j max /2. 
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Figure 20. Idealized Microburst: Simulation results of rain and hail at approximately 200s. hail (kg/m 3 ) 
using LeonardC scheme (top left); hail (kg/m 3 ) using MUSCL scheme (top right); rain (kg/m 3 ) using 
LeonardC scheme (bottom left); and hail (kg/m 3 ) using MUSCL scheme (bottom right), x-z plane is shown at 
y = jw/2. Only a portion of the domain is shown in the x direction. 
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Figure 21. Idealized Microburst: Simulation results of rain and hail at approximately 400s. hail (kg/m 3 ) 
using LeonardC scheme (top left); hail (kg/m 3 ) using MUSCL scheme (top right); rain (kg/m 3 ) using 
LeonardC scheme (bottom left); and hail (kg/m 3 ) using MUSCL scheme (bottom right), x-z plane is shown at 
y = jw/2. Only a portion of the domain is shown in the x direction. 
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Figure 22. Idealized Microburst: Simulation results at approximately 400s. hail (kg/m 3 ) using LeonardC 
scheme (top left); hail (kg/m 3 ) using the MUSCL scheme (top right); w-velocity (m/s) using LeonardC scheme 
(bottom left); and w-velocity (m/s) using the MUSCL scheme (bottom right). The x-y plane is shown at z = 
Zmax/2. Only a portion of the domain is shown in the x-y plane. 
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H. Tornadic Convection near Del City, Oklahoma 

The Del City, Oklahoma storm has been studied extensively in the past (e.g., Brandes 1981: Ray et al. 1981; 
Klemp et al. 1981; Proctor 1987). Beginning in the afternoon of May 20 th , 1977 and extending into the early hours 
of May 21 st , 1977, around sixteen storms associated with severe convection were observed in the Ft. Cobb, Edmond 
and Del City area in the state of Oklahoma. Two of the storms generated tornadoes (vorticity signatures of these 
tornadoes were observed by the Doppler radars deployed in the region) and supercells characteristics were observed 
during parts of their lifetimes. The estimated cost of damage from the tornados was around $600,000. In this 
section the results of two simulations using the Leonard scheme and the MUSCL-type scheme for the transport of 
microphysical scalars in TASS are presented for the Del City storm case. 

The computational domain was centered over Del City and bounded within [0:60km] x [0:60km] x [0:16.5km]. 
The mesh spacing was set to 600m in the horizontal and 275m in the vertical. The simulations were initialized (see 
Figure 23) homogeneously in the horizontal using a composite of two soundings from Ft. Sill (1500 CST) and 
Elmore City (1622 CST). The original soundings are shown in Figure 2 of Klemp et al. (1981). The composite 
sounding had a Lifted index of -8°C which may be lower than that of the actual environment. The wind fields for 
the initialization sounding were generated using observed wind tower data below 444m AGL and using a composite 
of winds from the Ft. Sill and Elmore City soundings for heights above 444m AGL. A spherical warm perturbation 
(AT = 1.5°C) with a radius of 10km was used to trigger the convection in the model simulation. The final time of 
the simulations was set to lOOmin. 

Figures 24-26 show the development of vertical velocity, hail and rain fields at approximately 24.17min, 35min, 
and 64.17min into the simulation respectively for the two schemes. It can be seen that although the macro-structure 
of the fields is quite similar even at 64.17min, differences start to appear as early as 24.17min into the simulation. 
This is a highly non-linear case in which the continuous feedback from the heat sinks and sources generated by the 
microphysical processes directly affects the development of the flow field. The differences in the magnitude and the 
distribution of the microphysical scalars due to different schemes starts altering the two flow fields accordingly - and 
these differences can be seen in the vertical velocity fields from the two simulations. 
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Figure 23. Del City, Oklahoma. Initialization sounding representative of conditions near Del City. 
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Figure 24. Del City, Oklahoma. Comparison of the MUSCL-type scheme with the Leonard scheme after 
1000 iterations (~24.17min into the simulation), x-z plane is shown in the figures. For hail and rain, the 
values represent the maximum along y direction, i.e., =ma x(q(i, j ,k)) . For the vertical velocity 

actual values for the x-z plane at y = y m *J2 are shown. The domain size is 60km x 60km in the horizontal and 
16.5km in the vertical. 
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Figure 25. Del City, Oklahoma. Comparison of the MUSCL-type scheme with the Leonard scheme after 
2000 iterations (~35min into the simulation), x-z plane is shown in the figures. For hail and rain, the values 
represent the maximum along j-direction, i.e., q(i,k) fi =ma x(q(i,j,k)). For the vertical velocity actual 

values for the x-z plane at y = j max /2 are shown. The domain size is 60km x 60km in the horizontal and 
16.5km in the vertical. 
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Figure 26. Del City, Oklahoma. Comparison of the MUSCL-type scheme with the Leonard scheme after 
5000 iterations (~64.17min into the simulation), x-z plane is shown in the figures. For hail and rain, the 
values represent the maximum along j-direction, i.e., q(i,k) f =ma x(q(i,j,k)). For the vertical velocity 

j=^,n 

actual values for the x-z plane at y = j max /2 are shown. The domain size is 60km x 60km in the horizontal and 
16.5km in the vertical. 
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IV. Summary 

Three upwind-biased schemes for the transport of microphysical scalars were analyzed and compared for 
implementation into the TASS model. All three schemes exhibit desirable properties however some differences 
were observed. The corrected Leonard scheme (LeonardC) has the better accuracy, and minimum diffusion and 
phase errors. The scheme however is not total variation diminishing and although undershoots are corrected, there 
is no mechanism to correct the overshoots in the solution. The MUSCL-type scheme and the high-resolution wave 
propagation method are both TVD and also exhibit minimal phase and diffusion errors. In terms of computational 
efficiency the MUSCL-type scheme is relatively faster compared to the LeonardC and the high resolution wave 
propagation method. Given the parallel constructs of the TASS model, the MUSCL-type scheme is also much easier 
to implement within the TASS model compared to the high-resolution wave propagation method. The MUSCL-type 
scheme was implemented into the TASS model and two simulations were conducted to evaluate the accuracy and 
computational efficiency of the scheme. The first case of an idealized microburst showed that the two schemes give 
comparable results for the microphysical transport. The differences due to the two schemes resulted in changes to 
the flow fields (the simulated velocity magnitudes were in general larger when MUSCL-type scheme was used). 
Finally the Del City storm was simulated. The results showed that although the simulated macro- structure of the 
storm was similar for the two schemes, differences at smaller scales appeared at later times which also altered the 
flow field simulated by the two schemes. In future, other higher-order TVD methods such as the weighted average 
flux (WAF) and the weighted essentially non-oscillatory (WENO) schemes (Toro 1999) can be evaluated for 
possible implementation into the TASS model. 
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